The Damping and Excitation of Galactic Warps by 
Dynamical Friction 

Robert W. Nelson and Scott Tremaine 
Canadian Institute for Theoretical Astrophysics 

ABSTRACT 

We investigate the dynamical interaction of galactic warps with the surrounding 
dark matter halo, using analytic perturbation theory. A precessing warp induces 
a density wake in the collisionless dark matter, which acts back on the original 
warp, transferring energy and angular momentum between the warp and halo 
(dynamical friction). In most cases dynamical friction damps the warp, but in unusual 
circumstances (such as a halo that rotates in the same direction as the precession of 
the warp, or a warp in the equatorial plane of an axisymmetric prolate halo) friction 
can excite a warp. The damping/excitation time is usually short compared to the 
Hubble time for realistic systems. Thus most warps cannot be primordial; they must 
be maintained by some ongoing excitation mechanism. 

1 Introduction 

The origin and maintenance of the 'integral sign' warps commonly seen in the outer parts of disk 
galaxies is a longstanding puzzle in galactic dynamics (for a recent review see Binney 1992). Many 
if not most disk galaxies — even isolated ones — appear to be warped, which implies that warps are 
either long-lived or repeatedly excited (Sanchez-Saavedra et al. 1990, Bosnia 1991). Following an 
initial suggestion by Lynden-Bell (1965), Hunter & Toomre (1969; hereafter HT) explored the 
possibility that warps are discrete modes of oscillation of a thin, self-gravitating disk. They found, 
however, that isolated disks with realistic density profiles actually support a continuum of bending 
modes; consequently, any packet of bending waves should disperse. Dekel & Shlosman (1983) 
and Toomre (1983) suggested a solution to this problem, which was analyzed in detail by Sparke 
& Casertano (1988; hereafter SC): realistic disks can support discrete bending modes if they 
are embedded in the static potential of a flattened dark matter halo — these modes are distorted 
versions of the neutral tilt mode of an isolated disk. Thus the gravitational influence of galaxy 
halos may play an important role in the maintenance of galactic warps; conversely, by studying the 
dynamics of warps in galaxies we may learn about the shape and dynamics of the unseen halos. 

Unfortunately, galaxy halos are not simply static potentials: the collisionless halo material — 
whatever its nature — must respond to the gravitational field of an embedded disk. A warped, 
precessing disk induces a density wake in the dark matter halo and the gravitational field from the 
wake will then act back on the disk; in this manner energy and angular momentum are transferred 
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between disk and halo in a manner analogous to the Landau damping of waves in a collisionless 
plasma or dynamical friction on a point mass travelling through a stellar system. In general, 
then, any bending mode of a galaxy disk is not precisely neutral; the question is whether the 
dynamical coupling between the halo and the disk is strong enough to damp — or excite, as we 
shall see — bending modes of a galaxy on timescales less than the age of the galaxy itself 1 . 

This issue was first investigated by Bertin & Mark (1980), who approximated the warp as 
a tightly-wound (WKB) bending wave, embedded in an infinite, homogeneous, and isotropic 
'halo' having a Maxwell-Boltzmann velocity distribution. Their analysis predicted that the 
amplitude of the bending wave would grow inside the corotation radius — and on a short timescale. 
Unfortunately, this interesting conclusion is suspect, because Bertin and Mark neglected the 
unperturbed halo's contribution to the vertical restoring force (see Section 3.1). 

In this paper we examine the role of a responsive halo in damping or exciting bending 
disturbances in galactic disks. We extend Bertin and Mark's seminal analysis to include warps and 
halos with more realistic shapes, as well as halos with anisotropic velocity distribution and net 
rotation. In particular, if the orbits of halo particles — we shall refer to these as 'stars' although 
their nature is unknown — are quasiperiodic, energy and angular momentum exchange between disk 
and halo occurs via stars whose orbits are nearly resonant with the frequency of the time- varying 
disk potential (Tremaine & Weinberg 1984). In this situation, the direction of energy and angular 
momentum flow can be from disk to halo or vice versa, so that dynamical friction with the halo 
may either damp or excite warps, depending on the details of the halo structure. 

Our analysis of this phenomenon is more elaborate than any in the literature so far, but still 
has a number of limitations. We do not consider the self-gravity of the halo response. We calculate 
the unperturbed orbits of stars in the disk plus halo system by assuming that the potential is 
spherical, and assume that the halo is axisymmetric for the purpose of computing the precession 
rate and energy of the bending disturbance; whereas cosmological N-body simulations of the 
collapse of protogalactic density peaks typically result in triaxial halos. We use linear perturbation 
theory, whose validity is suspect since warps often have substantial amplitudes. 

Despite these shortcomings, we believe that the principal results of our analysis are robust. 
In particular, we find that dynamical friction from the halo can either damp or excite bending 
disturbances in a disk, although damping is the most common outcome for plausible disk/halo 
configurations. The characteristic growth or decay time is generally less than a Hubble time, 
implying that dynamical friction plays a major role in determining the present properties of galaxy 
warps. 

In Section 2 we review the dynamics of bending disturbances in a thin disk embedded in 
a static halo potential. In Section 3 we relax this assumption and consider the coupling of the 



lc One key question about this kind of proposal is how long the central disk can thereafter remain tilted despite 
dynamical friction against the doubtless flabby halo' (Toomre 1983). 
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bending perturbation to the dynamics of the halo stars. In particular, we give general formulae 
for the damping/excitation rate of bending disturbances from dynamical friction with the halo. In 
Section 4 we present numerical calculations of damping/excitation times for reasonably realistic 
halo and disk potentials and warp shapes. Finally, in Section 5 we discuss the implications of our 
results. 



2 Warp Dynamics 



We review the behavior of linear bending waves of a zero-thickness, axisymmetric, self-gravitating 
disk embedded in a static halo potential, delaying until Section 3 consideration of the coupling of 
bending waves to a responsive halo. Many, but not all, of the results we describe are presented in 
HT, Sparke (1984), and SC. 



2.1 The equation of motion 



We work in cylindrical coordinates (R, (f>, z) and take the unperturbed disk to lie in the z = 
plane, so that position in the disk is described by the vector R = (R,(f>). We consider a razor-thin 
axisymmetric disk of surface density fJ,(R), rotating at angular speed 0(_R) > 0. The equation of 
motion in the inertial frame for small vertical displacements, Z(R,(f>,t), is (HT) 

T£ s (£ +a ik)' z = F > +F » + F ~" (1) 

where the three terms on the right side are the vertical components of the force per unit mass 
arising from the gravitational field of the static halo (Fh) and the disk (F sg ), and from any other 
sources (F ext ). We assume that the halo is axisymmetric and that its density is smooth near the 
disk, so that the halo potential can be expanded in a Taylor series 

02 TT 

U h (R,z) = U h (R,0) + \vl(R)z 2 , where v 2 h (R) = ~j^(R, z = 0) (2) 

is the square of the frequency of small vertical oscillations in the field of the halo. Thus 

F h (K,t) = -4(R)Z(K,t). (3) 

It is sometimes useful to decompose the radial force —£l 2 R into its contributions from the halo 
and the disk, 

tf(R) = tf h (R) + tf d (R), (4) 



where 

0?(Jf\= „„ JliW _ 

R OR n dK ' RdRj IR-R 
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The vertical restoring force from a slightly distorted disk is given by HT, 

.\Z(K,t)-Z(K',t) 



F sg (K,t) = -Gj dK'^R') 1 v |R _ R/ | 3 ' • (6) 

We may write the disturbance Z as a superposition of eigenfunctions of the form 

Z(R, t) = Re [h{K)e 1 ^-^ } , (7) 

where we assume m > without loss of generality. If u is real and m ^ 0, the pattern speed 
fi p = uj/m is the angular speed of the rotating frame in which the disturbance is stationary. 

When external forces are absent {F ext = 0), the dynamical equation (1) for a single 
eigenfunction reduces to 

\[u - mtt(R)] 2 - v 2 h {R)} h(R) = G dR'R'fi(R')[h(R)H(R, R 1 ) - h(R')I m (R, R')] (8) 

where 2 

[ ' ' Jo (R 2 + R> 2 -2RRi cos t/>) 3 / 2 ' Jm ^'^ i (jR 2 + E /2 _ 2RRI cos ^ )3/ 
these are Laplace coefficients and can be evaluated in terms of elliptic integrals. 

The component of the torque on the disk parallel to the z = plane may be written 
r = T x e x + T y e y , where 

Tx + ir y = [ dKfi(R)[F z (y - ix) - z(F y - iF x )}; (10) 



2 , (9) 



here F = (F x , F y , F z ) is the force per unit mass on the disk at R, composed of a radial component 



-£l 2 R and a vertical component — v\Z + F sg . Thus 



r x + ir y = -i j dRii(R)Re^[(tt 2 - vl)Z(R,t) + F sg ]. (11) 
Since the torque from the disk on itself must be zero we have (Kuijken 1991) 

= J dRii(R)Re^[tt 2 d Z(R,t) + Fsg]; (12) 
(a result which can also be derived by manipulating eqs. 5 and 6) and 

Tx + ir y = -i J dRfj,(R)Re i ' t '(nl - vl)Z(R,t). (13) 



2 The functions H and I m diverge as R' — > R but the integral in equation (8) is nevertheless well-defined. 
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These equations degenerate to = when ra / 1 but when m = 1 we may use equations (1) and 
(7) to reduce equation (12) to the form 

/•oo 

J dRR 2 ii(R)[tt 2 p - 2tttt p + (tt 2 h - v 2 h )]h{R) = 0. (14) 

This constraint provides a quadratic equation for the pattern speed O p in terms of the eigenfunction 
h(R); its unique feature is that the self-gravity of the disk does not have to be evaluated explicitly. 
In general, one of the roots of equation (14) is spurious and we shall assume that the correct root 
is the one with the smallest absolute value. The motivation for this choice is that we focus on the 
distorted tilt mode, which has a small pattern speed. 

Numerical solutions of the eigenvalue equation (8) for various disk and halo density profiles 
are given by HT, Sparke (1984) and SC. The following special cases offer useful insight: 



1. In the absence of disk self-gravity, equation (8) reduces to 

uj = mVL(R) ± Vh(R), (15) 

as does (14) when m = 1. In general, the combination of terms on the right side will 
depend on radius, so the eigenfunctions are singular (van Kampen modes), and any spatially 
extended initial disturbance will damp by phase mixing. 

2. If the halo is spherical, the disk is neutrally stable to tilting. To see this from equation (8) 
we use equations (5) and (9) to write 

nl(R) = — dR'R'n(R')[RH(R,R')- R' h(R, R% (16) 
R Jo 

it follows that when m = 1 and the halo is spherical (yh = fi/J a solution of (8) is u = 0, 
h(R) ex R (the 'tilt mode'). 

3. For a nearly spherical halo, the tilt mode is distorted and acquires a small but non-zero 
pattern speed. In equation (14) we may therefore drop the term proportional to 0^; 
moreover, in a first approximation we can replace the eigenfunction h(R) by the tilt mode 
h(R) ex R. Thus we obtain (SC) 

J ° o (n 2 -p 2 h )fi(R)R 3 dR 
p ~ 2l^n(R)^R)R 3 dR ' [ ] 

note that the angular frequency in the denominator includes the contribution from both 
disk and halo, while the frequencies in the numerator represent only the halo contribution. 
Equation (17) simply says that the precession rate of the warped disk is the torque from the 
halo divided by the off-axis angular momentum (Kuijken 1991); we shall refer to it as the 
Sparke-Casertano formula for the pattern speed. 
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The distorted tilt mode is closely related to the Laplace invariable plane 3 : each annulus of 
the disk is in the Laplace plane determined by the halo and the other annuli. Thus we shall refer 
to this mode as the 'Laplace mode' (SC use the term 'modified tilt mode'). 



2.2 Energy and angular momentum 



To determine the energy E w associated with a bending disturbance Z(R, /), we assume that the 
disturbance has been generated by an external potential U ex t(r,t), where r = (R,(f>,z) = (R, z). If 
the disk has density p(r,t) and velocity v(r,/) then the rate at which the external potential does 
work on the disk is 

= J drpv ■ F ext = - J drpv ■ VU ext . (18) 

We may write v = (0,fii2, DZ / Dt); the radial and azimuthal velocities have been set to their 
unperturbed values since we are working to first order in the small quantity Z(R, /) and to this 
order the horizontal and vertical motions are decoupled. Thus 

dE 

' drp [9,U ext ,4, + (Z, t +SlZ, <t) )U ext , z ] , 



dt 

dr [Slp^Uezt - pZ, t U ext , z + Slp, z Z,4> U ext ] , (19) 



where the second line follows from integration by parts. The density of our razor-thin disks may 
be written p(r,t) = p(R)S[z — Z(R, /)]. Thus 

p,0 = -pS'(z - Z)Z,j, = -p, z Z,4, . (20) 

Substituting this result in equation (19) we find that the first and third terms cancel, so that 

= - J drpZ, t U ext)Z = j dKpF ext Z, t (21) 

where as usual F ext is the ^-component of the external force (note that the expression is not 
/ dR,pF ex t(D Z / Dt) , as one might expect from considering only vertical motions). 

We may eliminate the external force F ext from equation (21) using equations (1), (3) and (6). 
Then using integration by parts with respect to (f>, we can convert the right side of equation (21) 



3 Consider a flattened planet that is accompanied by a massive satellite on an inclined orbit. A test particle 
orbiting the planet in the Laplace plane will remain in that plane (except for short-period oscillations) despite the 
perturbations from the planet's equatorial bulge and the massive satellite. The normal to the Laplace plane lies in the 
plane determined by the spin and orbital angular momentum vectors of the planet and satellite, while the orientation 
of the plane is fixed so that the angular momentum vectors of the planet, satellite and test particle all precess at the 
same rate. 
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into the total time derivative of a quantity which must equal the energy E w to within a constant. 
Setting this constant to zero when Z = 0, we obtain 

E w = \j dR^R)[Z 2 W h Z 2 - n 2 Z 2 ] + | / rfRrfRV(igK^) [Z(R |^:^ /,0]2 - (22) 

In the absence of an external force, the energy (22) is conserved 4 . 

Note that the energy of a disturbance is not guaranteed to be positive, since the first 
integrand can be negative. HT examined the sign of the energy for azimuthal wavenumbers m = 
(axisymmetric) and m = 1 (tilt or warp). When m = Z^ = so equation (22) implies that the 
energy is positive (or zero, in the trivial case = 0, Z =constant, corresponding to the translation 
of an isolated disk). Since unstable modes must have zero energy, all m = modes are stable. 

In the case m = 1 it is useful to decompose the angular speed into contributions from the disk 
and halo (eqs. 4, 5). HT then show that 



■\ J dRfi(R)tt 2 d (R)Z,l+- J dRdR'n(R)n(R 



|R- R'l 3 



= j J dRdR'd<f>n(R)n(R')[RZ(R', <f>, t) - R'Z(R, <f>, t)] 2 h(R, R'), (24) 
where I m is defined in equation (9). Thus the energy of an m = 1 disturbance may be written 

E w = \ J dRii(R)[Z, 2 +(ul-n 2 h )Z 2 ] + ( ^ j dRdR'd4>^{R)^{R')[RZ{R',4>,t)-R'Z{R,4>,t)] 2 I 1 {R,R'). 

(25) 

If there is no halo (yh = 0^ = 0) the energy is positive or zero; zero energy corresponds to the 
disturbance Z = ARcos(cj) — 4>q) where A and (f>o are constants (physically, an isolated disk is 
neutrally stable to a uniform tilt). Since unstable modes must have zero energy, all m = 1 modes 
of an isolated disk are stable, a conclusion due to HT. 

When the disk is embedded in a halo, m = 1 disturbances can have positive or negative 
energy. We shall call an axisymmetric halo 'vertically stiff' if Vh > 0^ at all radii (which is 
generally the case if the disk lies in the equatorial plane of an oblate mass distribution), and 
'vertically soft' if Vh < 0^ (which occurs for prolate distributions). Equation (25) shows that the 
total energy of an m = 1 disturbance is positive — so that m = 1 modes are stable — whenever the 
halo is vertically stiff. 



4 This result for the energy is perhaps somewhat unexpected: in the absence of self-gravity, the most natural way 
to derive an energy is to treat the disk particles as independent vertical harmonic oscillators. This procedure yields 

E' w = \ J dRn(R)[(DZ/Dt) 2 + vlZ 2 ] = \ J dK^(R)(Z, 2 t +v 2 h Z 2 + U 2 Z,% +2VZ, t Z,* ), (23) 

which is not the same as the first term of (22). Both E w and E' w are conserved if F sg = F ex t = 0, but equation (22) 
is the only energy-like quantity conserved when self-gravity is present (except in the special case of uniform rotation, 
which is described at the end of this section). 
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The expression for the energy (22) can be rewritten in a simpler form. We operate on equation 
^1) with f dR,p(R); setting F ext = we obtain 

J dRfi(R)(ZZ, tt +2nZZ, t<t> + n 2 zz,^W h z 2 ) 

, Z(R,t)[Z(R,t) - Z(R',t)] 



G / dRdR! n(R)n(R' 



|R- R'| 3 

= -f/^RVW^ [Z(R ^I^° ]2 , (26) 

where the last expression is obtained from the second by interchanging primed and unprimed 
variables and averaging with the original expression. Using this result to eliminate the last term 
from equation (22) we find 

E w = \ j dRp(R)[Z 2 -ZZ, tt -2QZZ, t<t> -tfZ,\ -tfZZ, H }. (27) 

The terms proportional to fi 2 cancel upon integrating by parts with respect to (f>. Thus 

E w = \ j dRp(R)[Z 2 -ZZ, tt -2QZZ, t<t> ]. (28) 

For a disturbance of the form (7) with m/0 and real frequency u = mO p , we have 

E w = -Trm 2 tt p J dRR[ji(R)\h(R)\ 2 (Sl - fi p ). (29) 

Thus disturbances with zero pattern speed have zero energy; if the pattern speed is negative 
(retrograde precession) the energy is always positive; if the pattern speed is small and positive the 
energy is negative. 

We can also work out the ^-component of angular momentum associated with a disturbance 
excited by an external potential: 

dJ i = -J dr P u °*'* 

= Jdr P ,,U ext 

= -jdr P , z Z,,U ext 

= J drpUext^Z,^, (30) 

where the third line follows from equation (20). We eliminate U ex t )Z = ~F ex t using equations (1) 
and (3) and derive an expression for the angular momentum by a similar procedure to the one 
used to derive equation (22). In this derivation we use the results 

Zrft = (Z^ Z, t ),t —\{Z 2 ),0 (31) 
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and 

. Z^{Z - Z') 



J dRfj,(R)Z, <t> F sg = -G J dRdR! n(R)ii(R' 



|R- R'l 3 



G f^^,.,„,.,„,d(Z-Z 



r\2i 



dRdR'n(R)n(R 



2 J rK Jrx ' |R- R'| 3 

~2 



J dRdK'fi(R)fi(R')(Z - Z') 2 -^\K - R'| -3 , (32) 



which is zero, because the integrand to the left of d/d(f> is symmetric under interchange of primed 
and unprimed variables, while 9|R — 'R'\~ 3 /d(f> is antisymmetric. 



In this way we find 



J w = - I dRii(R)(Z, t Z^ ). (33) 



Thus, in the absence of external forces, there are two independent conserved quantities that 
are quadratic in the disturbance strength. These quantities correspond to the energy E w and 
^-component of angular momentum J w . For a disturbance of the form (7) with real frequency 
u = mfip, we have 

J w = -7TTO 2 / dRRfi(R)\h(R)\\n - ftp). (34) 



Comparison of equation (34) with equation (29) shows that 

E w = Q,pJ w . (35) 

We close this section by considering the special case of disks with constant angular speed, 
0(_R) = fio =const. Combining equations (22) and (33) we obtain the conserved quantity 

E w - tt J w = \ j d^(R)[(Z, t +tt Z A f + v 2 h Z 2 ] + || dKdK'^RMR') [Z(R |^ ~ ^ )] ' . 

(36) 

Since this quantity is non-negative, all bending disturbances in a uniformly rotating disk are 
stable. 



2.3 The WKB approximation 



The equation of motion (1) for bending waves can be solved analytically in the WKB or 
short-wavelength limit. This solution offers useful insights, even though the WKB approximation 
is poor for the bending waves actually seen in galaxies. 



In the WKB limit the radial wavelength A is much shorter than any relevant radial scale 
length, i.e., 

din fJ,(R) 



A < R, A < 



dR 



(37) 
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but long compared to the wave amplitude, 

A>|Z|. (38) 

We write the vertical disturbance as 

Z(K,t) = Re !A(R,t)e t I RK{R ^ R e t{m ^'- ujt A , (39) 



where K(R) = 2ir/\(R) is the local wavenumber and A(R,t) is a slowly varying amplitude. 

The vertical gravitational force from the disturbed disk (eq. 6) can be evaluated using local 
Cartesian coordinates, R' = R + xe x + ye y , where at R the a; — axis points radially outward and 
the y— axis points azimuthally. Because the radial wavelength is short, we neglect the radial 
variation in the surface density fJ,(R') and the azimuthal variation in Z, and write 

Z(R',t)= Re[Aoe^ Kx -^], Z(R,t) = Re^oe"^], (40) 

where 

A = A(R,t)e l f RK ^ dR e tm ^ R \ (41) 

Then equation (6) yields 

F. g (R,t) = -CMR) Re { Aoe-*" j dx dy ^'^^ } 

/ \ i i / x f°° dv /» , (1-cosm) 

= -4C M RU\z(R. t ) I j^-^l 

= -2TrGn(R)\K\Z(R,t). (42) 

The WKB dispersion relation follows from the dynamical equation (1), 

D(u, K) = [u- mtt(R)] 2 - 2nGfi(R)\K\ - v 2 h {R) = 0. (43) 

At the corotation radius, the angular speed of the disk equals the pattern speed, 0(_R) = fi p = uj/m. 
At the vertical resonance radii K(R) = 0, that is, u = raO(_R) ± fh{R) (cf. eq. 15). Because 
the dispersion relation can only be satisfied when m 2 (i7 — fi p ) 2 > there is a forbidden region 
between the two vertical resonance radii and including corotation, in which WKB bending waves 
cannot propagate. 

The radial group velocity for WKB bending waves is (Toomre 1969, Whitham 1974) 

dcj irGfi(R) sgn(ii') 
Cg = ~dK = m[tt(R)-tt p ] ' 



(44) 



Thus the direction of propagation depends on whether the waves are leading (K < 0) or trailing 
(K > 0), and on whether they are inside corotation (0 > O p , assuming that the angular speed 
decreases outward) or outside. 
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The equation of motion for vertical disturbances (1) can be derived from a Lagrangian. The 
dynamics in the WKB limit are therefore described by an averaged Lagrangian that is obtained 
by averaging over the short-wavelength oscillations (Whitham 1974), 

L(u, K, A) = IttRi2\A\ 2 D(uj, K), (45) 

where |A| is the wave amplitude and D(uj,K) is the dispersion relation (43). The energy density 
of the wave per unit radius is 

OL 

e(R) = uj— = -Trm 2 Ri2(R)n p (tt - tt p )\A\ 2 . (46) 

The energy density is positive if the pattern speed is prograde (0 p > 0) and outside corotation 
(0 < fip), or if the pattern speed is retrograde. The total energy of the disturbance is 

E w = J e(R)dR, (47) 

which is the same as the energy derived in the previous section (eq. 29) — in this case it happens 
that the WKB approximation is exact. 

The density of z— angular momentum per unit radius is (Bertin & Mark 1980) 

, , dL e(R) , s 

« = ^ = V' (48) 

which is exactly consistent with equation (34). 

3 Dynamical Friction Between a Warped Disk and a Spherical Halo 



So far we have treated the halo potential as a static background field, which influences the 
disk only through its contribution to the vertical and radial restoring force. In fact, the halo 
is composed of stars whose orbits respond to the time- varying gravitational field of the warped 
disk. The gravitational attraction between the coherent response of the halo and the warped 
disk transfers energy and angular momentum between the two components, thereby damping — or 
possibly exciting — the warp (dynamical friction). 

To describe this response, we write the halo potential as Uh(r) = U®(r) + £7^(r,/), where the 
two terms denote the static, unperturbed halo potential and the potential arising from its response 
to the warp. The corresponding densities are /0°(r) and p\(r,i). For consistency with Section 
2.2, we define the energy of the bending disturbance in the disk, E w , to include the kinetic and 
self-gravitational energies of the disk material plus its potential energy in the static field U®, but 
not its potential energy in the response field U^; the forces arising from the latter field are treated 
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as an external force F ext = — VC^. Thus 

dE,„ 



dt j drp d v d ■ F ext = J drV ■ (p d v d )Ul 

. d Pd TT i _ r , f j^dpd p{ 



dv^Ul = G / drdr' , 

dt h J dt \r - r'\ 

x dU d f x dU w 
dr Ph — = - J dr Ph — . (49) 

where we have integrated by parts and used the continuity equation for the disk mass; the subscript 
'd' denotes disk, and in the last line we have replaced the disk potential U d by the warp potential 
U w since the potential of the unwarped disk is stationary. To first order in the perturbation 
strength we may write p\ = — V • (p®Ar), where Ar(r,i) is the Lagrangian displacement of the 
halo stars caused by the perturbation. The last line of equation (49) may now be written 

dEin f * n . -dU v 



dt 



jdrp°Ar.V^. (50) 



There are two equivalent ways to determine the dynamical friction exerted by the halo on 
the disk: (i) calculate the density and potential perturbation induced in the halo by the disk and 
then integrate the resulting gravitational force from the halo over the disk mass distribution — this 
is analogous to the standard discussion of Landau damping in a collisionless plasma (c/. Ichimaru 
1973) and corresponds to the first line in equation (49); (ii) calculate the rate of change of energy 
of each halo star due to the warp, and integrate these energy changes over the halo (which 
corresponds to eq. 50). We use the second approach in this paper. 

Our calculations are based on the action-angle formalism developed in papers by Lynden-Bell 
& Kalnajs (1972), Tremaine & Weinberg (1984), and Weinberg (1986), which may be consulted 
for details. Our analysis employs three (not very good) approximations: (i) the warp is small, so 
that we can use perturbation theory in which the small parameter is the amplitude of the warp; 
in this approximation both Ar and dll d /dt are proportional to the warp amplitude so dE w /dt 
is second-order in the amplitude (this is the same order as the energy of the warp computed in 
Section 2.2; thus the damping rate E w /E w is independent of amplitude); (ii) the unperturbed 
potential (halo plus flat disk) is spherical; this is certainly incorrect — and moreover inconsistent, 
since we assume that the halo is axisymmetric but non-spherical when estimating the pattern 
speed and energy of bending modes — but orbits in spherical potentials are easier to analyze than 
orbits in flattened potentials and the dynamical friction force should not depend strongly on 
the shape of the halo orbits; (iii) we neglect the self-gravity of the halo response; again, this is 
certainly incorrect but unlikely to cause more than a factor of two or so error in the results. 

In the unperturbed state (no warp) the motion of the halo stars is described by an autonomous 
Hamiltonian Hq. The phase-space coordinates of a star are specified using the action-angle 
variables of this Hamiltonian (I,w); in these coordinates Hq = Ho(I) so that Hamilton's equations 
become 

I = "^ = 0, w=^Ufi(I), (51) 
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with solutions 

= Io = constant w (0 = w o(0 = fi(Io)^ + constant. (52) 

The density of halo stars in phase space is described by a distribution function (hereafter df) F 
and Jeans' theorem implies that F = F(l). 

When a warp is present the Hamiltonian becomes 

H = H + U w . (53) 

Hamilton's equations become 

1=-^—, w = n(I)+— — . (54) 
The perturbing potential can be expanded in an action- angle Fourier series (cf. Appendix A), 

{oo oo I 

£ £ Vwsiiy^-"!*]. (55) 
/ 3 =0/i,/ 2 = -oo J 

The first-order corrections to the orbit (52) caused by the perturbation are 

AI=-|^, Aw = | (56) 
ow ol 



where 

rt f oo oo e i(l-W-wjt) 1 

X = j_ oo U w dt = RelYE -<») ]> (57) 



and these quantities are evaluated along the unperturbed orbit (Iq,wo). We have assumed that 
uj\ has a small positive imaginary part, tj, as if the potential were turned on slowly in the distant 
past. 

The integrand of equation (50) may be written 

Ar - V ^r = S« I — «>AI + £^(l.,wo,0Aw. < 58 » 

We evaluate this expression using equations (55)-(57). Since the df depends only on the actions 
we can average the result over the unperturbed angles wo- As described in Tremaine & Weinberg 
(1984), the terms in this expression involve combinations of periodic exponentials, most of which 
vanish upon averaging over initial phases 5 . Combining the surviving terms we find 

Ar-V^\ =f f 1- * W)|» ^l- 1 -")^ \ (59) 

a* /w hu,t-J 1( jl 2(|i- n-^ + ^)/' (59j 



5 Or are oscillatory in time, in the case k = 0. 
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where lo\ is now taken to be real as the contribution from the small imaginary part t] is written 
out explicitly. We now use the identity 

lim — = irS(x). (60) 

Integrating over the df of halo stars finally gives the rate of change of the disk energy 6 (cf. eq. 50), 
E d = -JdIdwF(I)Ar.V^ = -(27r) 3 JdIF(I)^Ar.V^ 

= 4vr 4 £ £ ^1(1.0)1.^1^(1)1^(1.0-^) 



/ 3 =0 li ,l2 = — CO 

CO 



" 47r4 E E J ^.1(1.0)^(1)1^(1)1^(1.0-^), (61) 

where dS is the area element on the boundary of action space and we have set / = 0. For 
most dfs the last (surface) term vanishes, as the following argument shows: Let dS be an 
area element on the boundary of action space. The flux of stars through this element with 
angles in the interval dw is dwdS ■ I-F(I); using equations (54) and (55) this flux equals 
— J2l dwdS ■ l_F(I)Re{i$i(I) exp[i(l • w — The terms in this sum are linearly independent 

functions of w. Since stars cannot leave action space, each term in the sum must be zero. Thus 
dS ■ l_F(I)$i(I) must vanish, which implies that the surface term in (61) must vanish. (This 
argument can fail if the df diverges on the boundary of action space.) If the surface term vanishes, 
equation (61) simplifies to 

Ed = ^ 4 jt £ /dI(l-n)l-^Pl*i(I)| 2 «(l-n-^i). (62) 

Only resonant orbits with 1 • O(I) = u\ contribute to the energy transfer; a particular resonant 
triplet 1 = (h,l2,h) contributes to a net damping or growth of the perturbation depending on the 
components of O(I) and the gradient of F(l) along the resonant triplet. 

Finally, we note that the formalism we have developed is related to Goodman's (1988) stability 
criterion for galaxies, which is based on the sign of a quadratic functional of a trial perturbation. 
Consider a trial perturbation U w with time dependence with rj real (not e - im ^ t v t + r i t as we have 
assumed so far). Then Goodman shows that a sufficient criterion for instability is 



r\G f , , / , x , ,JZ(K,t) - Z(K't)r , s 

E w > J yJ dRdR'fi(R)fi(R') 1 { ft_ R ,| 3 ^ , (63) 



6 Our result is consistent with equation (65) of Tremaine & Weinberg (1984), who assume k ^ and a single pattern 
speed, so that Lo-y = h^lp and the torque is related to the rate of energy change by E = fi p r z (Jacobi's integral). 
We have assumed that the evolution of the warp is sufficiently rapid that we are in the 'fast' limit of Tremaine and 
Weinberg, that is, that the damping rate is large compared to the perturbation strength; this assumption is always 
justified if the warp amplitude is sufficiently small, because the damping time is independent of amplitude. 
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where E w is defined by equation (49) and evaluated using equation (59) — compare Goodman's 
equation (24). 

3.1 Halos with isotropic distribution functions 

An important special case is an unperturbed halo df that depends on energy alone, 
F = F(E) = F[H (I)]. Then equations (51) and (62) imply 




In most model galaxies with F = F(E) the halo df is a decreasing function of energy, (IF I (IE < 0; 
following Goodman (1988), we call these 'iddf halos' (for 'isotropic, decreasing distribution 
function'). In iddf halos, equation (64) implies that E,i < 0. In other words, in iddf halos, a 
warped disk always loses energy to the halo. Thus disk disturbances with positive energy damp, 
while disturbances with negative energy grow. 

Equation (22) shows that bending disturbances with a single azimuthal wavenumber m > 
have positive energy if i>h{R) > mCl(R) at all radii. Thus in iddf halos, dynamical friction damps 
all bending disturbances so long as 

v h (R) > mSl(R) at all R. (65) 

For m = 1 disturbances, equation (25) implies that the damping criterion (65) can be replaced by 
the stronger damping criterion 

v h (R) > Sl h (R) at all R, (66) 

which is the condition that the halo is vertically stiff. In other words, in vertically stiff iddf halos, 
all m = 1 bending disturbances damp. 

Bertin & Mark (1980) examined the evolution of short-wavelength bending waves of a 
thin disk embedded in a halo. The halo was assumed to have an isotropic Maxwellian velocity 
distribution, and the disk dynamics were analyzed in the WKB approximation. They concluded 
that dynamical friction from the halo would excite all bending waves propagating inside the 
corotation radius. Toomre (1983) has pointed out that the validity of this conclusion is suspect, 
because Bertin and Mark ignored the halo's contribution to the vertical frequency, effectively 
setting Vh = and thereby ensuring that the halo was vertically soft. The more general approach 
developed in this paper allows a rapid derivation of the correct answer: if, like Bertin and Mark, 
we assume an iddf halo, then bending waves are excited by dynamical friction from the halo if 
and only if they have negative energy. From equation (46) the energy density of a WKB bending 
wave is negative if and only if 

< ft p < Sl(R), (67) 
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i.e. for prograde waves inside corotation. 

Equation (67) is the condition for excitation of WKB bending waves for any iddf halo. The 
WKB dispersion relation (43) for m = 1 waves implies that |fi — p | > z/^, so a necessary condition 
for excitation is < fi, which is consistent with the damping criterion (65) 7 . These results imply 
that the excitation of WKB bending waves by the halo is considerably more difficult than Bertin 
k Mark (1980) assumed. 

We shall not pursue the analysis of WKB bending waves further, since the WKB approximation 
is invalid for most observed galactic warps. Instead we shall focus on numerical evaluation of 
damping and excitation rates for more realistic warps. 

3.2 A model for the halo 

In this section we describe the analytic model for the halo df that we shall use in our numerical 
estimates of damping rates. 

In a universe dominated by collisionless dark matter, galaxy halos form by the gravitational 
collapse of random peaks in the cosmological density field (c/. Warren et al. 1992 and references 
therein). Both analytic arguments and numerical simulations show that halos have power-law 
density profiles that are roughly consistent with the profiles p oc r~ 2 implied by observations 
of flat rotation curves of spiral galaxies. The simulated halos tend to be triaxial with minor to 
major axis ratios c/a ~ 0.5, although the uncertain effects of dissipation are likely to reduce the 
triaxiality (Katz k Gunn 1991; Dubinski 1994). The velocity ellipsoid in the halo is radially 
elongated: the ratio of the mean-square radial velocity to the mean-square tangential velocity is 
given by 2{v 2 ) / {v 2 ) ~ 1.4 (Dubinski 1992), whereas an isotropic dispersion tensor would have the 
value unity. Tidal torquing from neighboring structures gives the halo a net angular momentum; 
Dubinski finds that the ratio of the maximum rotation velocity to the central line-of-sight velocity 
dispersion is roughly correlated with the dimensionless spin parameter, f max /<7 c ~ 4A, where 
A = J\E\ 1 / 2 G~ 1 M- 5 / 2 has a median value of about 0.05 (White 1984; Barnes k Efstathiou 1987; 
Zurek et al. 1988). 

To simplify our calculations we ignore the intrinsic triaxiality of the halo and the contribution 
of the disk to the total potential: thus we assume a halo df, halo orbits, and a rotation curve 
for the disk that are consistent with a spherical potential provided entirely by the halo. This 
approximation is not strictly consistent, since we assume that the halo is axisymmetric but 
non-spherical when estimating the pattern speed of bending modes or calculating the mode energy 



r In the WKB approximation the damping criteria (65) and (66) are generally the same when m = 1, because 
there is no distinction between fi^ and fi. The reason is that the dispersion relation (43) can usually be satisfied 
in the WKB limit (\R ' R\ 1) only if the surface density fi is small, so the disk's contribution to the total rotation 
speed fi is negligible. 
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(see eq. 75 at the end of this subsection). The bias caused by the spherical approximation should 
lead us to underestimate the damping and excitation rates in the likely case that the halo is 
flattened towards the disk plane, since the halo density will be higher and the velocity of the halo 
stars will be lower than we have assumed. 

Our model for the halo df: (i) assumes that the halo potential generates a flat rotation curve 
with circular speed v c , 

U h (r) = v 2 ln ( -) , (68) 



which implies a halo density 



P(r) = = ^kr 2 (69) 



V 2 U{r) _ _v 2 
4ttG ~ AirGr 

(r s is a fiducial radius, which eventually scales out of the problem); (ii) allows for a range of 
velocity ansiotropy, expressed by the ratio 2{v 2 ) / {v 2 ); (iii) allows a non-zero mean rotation speed 
{v^), which is taken to be independent of radius. 

The df that we use is 

F(E,J,J Z ) = F J 2 ^e- E ^\l + aJ z /J), (70) 

where E is the energy, J is the angular momentum, and J z is the component of angular momentum 
normal to the unperturbed disk. The parameter a controls the halo rotation, while 7 controls the 
velocity anisotropy; we must have |a| < 1 so that the df is non-negative, and 7 > so that the 
spatial density is finite 8 . In order that this df generates the density (69) we must have 

* 2 - v l r 1 r< +1/2 (7Y) 

where J s = r s v c . In the special case a = 0, 7 = 1, the df and potential reduce to those of the 
singular isothermal sphere (the normalization is the same as Weinberg 1986, eq. 24). In turn, F is 
a special case of the larger classes of scale-free dfs discussed by Gerhard (1991) and Evans (1994). 

The mean-square velocities in the radial and tangential directions are 

{v 2 ) = a 2 , {v 2 ) = 2 7 cT 2 , (72) 

so the ratio of the principal axes of the velocity ellipsoid is given by 

ff-il (73) 



8 Galaxies with predominantly radial orbits are likely to be unstable; in fact Palmer & Papaloizou (1987) have 
shown that all galaxies with distribution functions that are unbounded as J — > are unstable. However their proof 
assumes that the resonant combination of frequencies Qi — 2Q2 oc J as J — > 0; while this is generally true for galaxies 
with a flat central core, in the singular isothermal sphere Qi — 2Q2 oc 1/| In J| as J — > 0, and the Palmer-Papaloizou 
proof does not apply. In any case, it is not crucial for our purposes to work with stable halo models, because the 
self-gravity of the halo response — which is central to the instability — is neglected in our calculations of dynamical 
friction. 
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a halo with nearly radial orbits corresponds to 7 — ► while one with nearly circular orbits has 
7 — ► 00. The mean rotation speed is 

^)_^_r(7+i) 

Vc -27 1 / 2 r( 7 ) sm "- [74) 

The cosmological simulations described above suggest that typical parameter values are 7 ~ 0.7 
(from the velocity anisotropy) and a ~ 0.4 (we have identified f max with {v^)(6 = ^tt) and 3a 2 
with (v 2 ) + (v 2 )). 

Finally, when estimating the pattern speed or energy of a bending mode, we generalize the 
expression for the halo potential (68) to an axisymmetric halo: 

U h (R,z)= lv 2 c ln[R 2 + z 2 /ql) , (75) 

where g$ is the axial ratio of the equipotential surfaces and e$ = 1 — g$ is their ellipticity (e$ > 
for oblate halos and < for prolate halos). In this case the vertical and azimuthal frequencies in 
the field of the halo are (eqs. 2 and 5) 

Vh = ^R> = (76) 

3.3 A model for the warp 

To calculate the damping rate in a realistic warped galaxy we need to know the shape of the warp 
h(R) (eq. 7) and its associated gravitational potential U w (r). Rather than solving equation (8) to 
obtain an exact eigenfunction, we have chosen a simple analytic parametrization of the warp. We 
assume that the surface density of the disk is exponential, 

fi(R) = ^oe- R/Rd , (77) 
so that the total disk mass M,i = 2iri2oR 2 . We set the parameter 



GM d 
R d v 2 c 



1.6, (78) 



in all of our simulations, which is roughly correct for our Galaxy. We write h(R) = Rdf(R/Rd) 
where Rj is the scale length of the disk. If the mode has real eigenfrequency then we can assume 
that h(R) is real (by the anti-spiral theorem of Lynden-Bell & Ostriker 1967); in this case the line 
of nodes is radial, which is approximately consistent with observations. 

Our parametrization has the form 

f{x) = Cx n e- Xx ; (79) 



where C is an arbitrary normalization (recall that the damping rate is independent of the 
amplitude). This parametrization permits us to model two distinct types of warps, (i) For n = 1 
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the inner part of the disk is tilted relative to the halo (as in the Laplace mode) and the outer disk 
is flat; in this case equation (79) gives a warp whose curvature in its outer parts is more gradual 
than observed warps (see Figure la), but this is not a major defect since most of the interaction 
with the halo comes from the tilted inner disk, which is not curved either in our model or in real 
galaxies, (ii) For n ^> 1 the inner disk is flat and the outer disk is tilted. In this case the shape of 
the warp can be made to closely resemble observed warps (except that at very large radii the disk 
becomes flat again — its maximum height is at x max = n/X — but if x max ^> 1 the surface density in 
the region where the disk becomes flat again will be negligible). 

An advantage of the functional form (79) is that the spherical harmonic expansion of the warp 
potential Ui n (r)Yi n (0 , (f>) can be evaluated in terms of incomplete gamma functions (see Appendix 

B). 
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Fig. la. — The solid line shows the shape of a typical warp used in our estimates of the warp 
damping/excitation rate, derived from eq. (79) with C = 1, n = 1, A = 0.2. The dotted and 
dashed lines indicate the first two non-zero components, U2\(r), U^iir), of the spherical harmonic 
expansion of the potential of the warped disk with m = 1; the scale for the potentials is shown on 
the right hand axis in units of GM^/Rd- 
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Fig. lb. — The pattern speed in units of v c /Rd (the angular speed of the disk at one scale length), 
as a function of the ellipticity of the halo potential. The pattern speed is computed using the warp 
shape in Figure la and eq. (14) — the solution of the quadratic equation with the smaller value of 
| Op | is shown. The dotted line corresponds to the Sparke-Casertano formula (eq. 17). 

Fig. lc. — Energy of the warp shown in (a) as a function of the ellipticity of the halo potential, 
computed from eqs. (29), (77) and (78). 
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Figure la shows an example of a warp shape parametrized by equation (79) with n = 1 and 
A = 0.2, along with the radial potential functions associated with the first two nonvanishing m = 1 
spherical harmonics, U2\(r), U^iir). Figure lb shows the pattern speed for this warp, calculated 
using equation (14), as a function of the halo ellipticity e$ (eq. 75). The dotted line shows the 
Sparke-Casertano formula for the pattern speed (eq. 17), which is a good approximation to the 
more accurate formula. Figure lc shows the energy of the warp (eq. 29) with C = 1. The figure 
illustrates that warps embedded in vertically stiff halos (e$ > 0) generally precess in a retrograde 
direction (O p < 0) and have positive energy, while warps embedded in vertically soft halos (e$ < 0) 
precess in a prograde direction (O p > 0) and have negative energy. 



4 Calculation of Damping and Excitation Rates 



We now calculate the rate of change of the warp energy, and the resulting damping or growth rate 
caused by dynamical friction from the halo. Weinberg (1986) has described similar calculations, 
based on a formula analogous to (62), for the rate of decay of a satellite on a circular orbit in a 
spherical potential. Thus we have simply summarized the principal formulae in Appendix A and 
refer the reader to Weinberg (1986) for details. 

The formula (62) for the rate of change of the warp energy is combined with equation (29) for 
the energy to obtain the characteristic evolution rate 

„ 1 (IE,,, 

r > implies excitation and T < implies damping. In our model the evolution rate 1 is a 
function of the velocity anisotropy of the halo (parametrized by 7), the rotation of the halo 
(parametrized by a), the warp pattern speed (written in dimensionless form as Q, p Rd/v c ), and the 
shape of the bending disturbance; it is proportional to the mass of the disk but independent of 
the amplitude of the disturbance (for small perturbations) since both the energy and the rate of 
energy change scale with the square of the amplitude. 

We evaluate the action-angle transform of the perturbing potential using equations (A1)-(A3) 
and (A23), along with equations (79) and (B9)-(B10) for the warp shape and potential. We 
treat the integrals in equations (A10), (All) and (A23) as the solutions to a coupled set of 
ordinary differential equations, which we integrate using an adaptive Runge-Kutta routine. We 
use Romberg integration to evaluate the integral in equation (A26), summing over the possible 
resonant triplets in the integrand itself, rather than evaluating an integral for each term separately. 
Since the number of resonant triplets grows oc / 3 , we must cut off the summation at a relatively 
small maximum harmonic / max . Experiments on a few points show that £(/ max ) = £ + C/^f x , with 
(3 ~ 1. In practice we evaluated £(/ max ) at / max = 2,4, and 6, fitting these values to a quadratic 
in l^ ax , and extrapolate to / max — ► 00. We believe that the damping/excitation rates we have 
obtained are accurate to within ±3%. 
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We have verified our results in the epicycle limit (1 — k <C 1) where the functional form for the 
potential transforms becomes analytic, and also compared our formulas with similar transforms 
calculated by Weinberg (1985). 
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Fig. 2. — The damping rate T = E w /E w as a function of the ellipticity of the halo potential. The 
warp shape is shown in Figure la. The halo df is described by eq. (70) with 7=1 (isotropic velocity 
distribution). Warps with T > grow, while those with T < damp. Damping or excitation occurs 
within a Hubble time if |T| > 0.0015 (assuming Rd/v c = 1.5 X 10 7 y). The small errorbar centered 
at e$ = corresponds to the range of damping rates where damping or excitation does not occur 
within a Hubble time. The dotted, solid, and dashed lines correspond respectively to halos with 
rotation parameters a = —1,0,1 described in the text; these values span the range of possible 
rotation rates. 

Figure 2 shows the dimensionless damping/excitation rate TRd/v c as a function of the 
halo ellipticity e$. The warp has the shape shown in Figure la and is embedded in a halo 
with an isotropic velocity distribution, 2(t^)/(ff) = 1. For parameters corresponding to those 
of our Galaxy, Rd/v c = 1.5 X 10 7 y, the damping/excitation time is less than 10 10 y when 
\T\R d /v c > 0.0015. 
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We plot results for three values of a which span the range of possible rotation rates consistent 
with our model (eq. 74). First consider the non-rotating halo (a = 0) with isotropic velocity 
distribution (an iddf halo in the notation of Section 3.1). As shown in that section, a warp 
always loses energy to an iddf halo. Thus warps with positive energy damp. However, if the 
halo is vertically soft (e$ < 0, which obtains if the disk lies in the equatorial plane of a prolate 
halo potential), the energy of the warp is negative (Figure lc). Thus, a warp embedded in a 
non-rotating, vertically soft iddf halo can grow in amplitude, while one embedded in a vertically 
stiff halo must damp. 

Halo rotation tends to increase the rate of change of warp energy E w if the halo rotates in 
the same sense as the pattern speed, and decreases E w if the halo rotates in the opposite sense. 
Thus, adding rotation tends to increase (or decrease) the warp amplitude if the sign of a£l p /E w 
is positive (negative). If the pattern speed is small, E w /Cl p is always negative (eq. 29); thus halo 
rotation contributes to warp growth only if the halo rotates in the opposite sense to the disk. 




Fig. 3. — Same as Figure 2, but with for a halo with predominantly radial orbits, 2(t^)/(ff) = 1.5. 
Note that the vertical range of the graph is a factor of ten larger than Figure 2. 

Figure 3 shows the damping/excitation rates for warps embedded in halos with predominantly 



radial orbits, 2(t^)/(ff) = 1.5. The damping rate is much larger than in the isotropic case; in fact 
so large (|T| ^ that our treatment of the warp as a slowly evolving normal mode is suspect. 




Fig. 4. — Same as Figure 2, but with for a halo with predominantly circular orbits, 2(t^)/(ff ) = 0.5. 
Note that the vertical range is only half as large as Figure 2. 

Figure 4 shows the damping/excitation rates for warps embedded in halos with predominantly 
circular orbits, 2(t^)/(ff) = 0.5. In this case energy flows from the halo to the warp, a consequence 
of the strong anisotropy in the phase space distribution function. Thus positive energy warps 
(e$ > 0) grow, while negative energy warps damp. 

Figure 5 shows the dimensionless damping rate as a function of the halo velocity anisotropy, 
for a warped disk with a retrograde pattern speed, R^p/vc = —0.5; this pattern speed corresponds 
to a halo whose equipotentials have ellipticity e$ = 0.5 (Figure lb). For halos with predominantly 
circular orbits the warp is excited, while in halos dominated by radial orbits it is damped. An 



interesting feature is that the damping rate diverges as [v 2 r ) / '(ff) — ► 2, corresponding to 7 




We show in Appendijx C that the apparent divergence is an artifact of the perturbation theory 
that we have used: according to equation (61) the rate of energy loss E w arising from the response 
of the halo to a perturbation potential $ is 0($ 2 ), but in fact the response is 0($ p ) where 
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Fig. 5. — The damping rate T as a function of the velocity anisotropy of the halo. The halo 
equipotentials are oblate, e$ = 0.5. The dotted, solid, and dashed lines correspond to halos with 
dimensionless rotation a = —1.0, 0.0, 1.0. Halos with predominantly circular orbits excite the warp, 
while halos with predominantly radial orbits damp the warp. 
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p = min(2, 27 + 1); thus for 7 < ^ the energy loss rate is of lower order than 0($ 2 ). Our analysis 
is not powerful enough to determine the damping rate for halos with 7 < ^; nevertheless it is clear 
that the damping is very strong when the orbits of the halo stars are predominantly radial. 

We have also examined the damping of disks that are only warped relative to the halo in their 
outer parts (eq. 79 with n = 5, A = 1). Although the damping/excitation times are substantially 
longer, significant damping still generally occurs in less than a Hubble time. 

5 Discussion and Conclusions 

The most striking feature shown by Figures 2-4 is the rapid evolution of the warps from dynamical 
friction. All of the models we examined had damping/excitation times less than a Hubble time 
when the halo ellipticity was larger than 0.2, and in general the damping/excitation times were 
much shorter than a Hubble time. For example, in a non-rotating halo with a plausible value of 
the velocity anisotropy, 2{v^)/{v^) = 1.5 (recall from Section 3.2 that cosmological simulations 
yield a typical value of 1.4 for this parameter), the damping time — T _1 = 10 8 y even when the 
ellipticity of the halo potential is as small as e$ = 0.2. 

Thus our numerical results show that dynamical friction between a warped galaxy disk and 
its surrounding dark matter halo usually leads to significant evolution of the warp over a Hubble 
time. The rate of evolution depends strongly on the velocity anisotropy and net rotation of the 
halo stars and the ellipticity of the halo potential. 

Dynamical friction can either damp or excite a warp: warps can have either positive or 
negative energy, and friction can either add energy to or remove energy from the warp. If the disk 
lies in the equatorial plane of an axisymmetric halo, then the energy of an m = 1 warp is generally 
positive if the halo is oblate or vertically stiff, and negative if the halo is prolate or vertically soft. 
A non-rotating iddf halo always removes energy from the warp. Halos with predominantly radial 
orbits remove energy more rapidly than halos with an isotropic velocity distribution, while in the 
unlikely case that the halo orbits are predominantly circular, dynamical friction can add energy to 
the warp. The effect of halo rotation is to add energy to the warp if the mean rotation rate is in 
the same direction as the precession of the warp, and to remove energy otherwise. 

Our assumption that the halo is axisymmetric is a major oversimplification, since cosmological 
simulations of galaxy formation typically produce halos that have highly tri axial shapes. The 
formation of gaseous disks in triaxial potentials has been investigated by analytic arguments 
and numerical simulations (Steiman-Cameron & Durisen 1984, Habe & Ikeuchi 1985, Pearce & 
Thomas 1991, Thomas et al. 1994). Disks can form in the principal plane normal to the major 
or the minor axis (orbits in the principal plane normal to the intermediate axis are unstable). 
As a first approximation, we may treat the halo as vertically stiff if the disk is normal to the 
minor axis and vertically soft if the disk is normal to the major axis. The principal plane in 
which accreted gas will settle depends on its initial angular momentum; simple analytic arguments 



based on perturbation theory imply that the probability that gas with randomly directed angular 
momentum will settle into the principal plane normal to the major axis is 2i cr i t /7r, where 



sin 2 i" crit = 



b~ 2 



— a 



-2 



(81) 



c 



-2 



— a 



Taking typical halo axis ratios of c/a = 0.45, b/a = 0.7 (Dubinski & Carlberg 1991, Warren et 
ai. 1992), we find that the probability that disks lie in the principal plane normal to the major 
axis is roughly 0.3; this is an upper limit since in practice the angular momentum vector of the 
infalling gas is strongly correlated with that of the halo material (e.g. Katz & Gunn 1991) For 
the minority of galaxies in which the disk forms in the principal plane normal to the halo major 
axis, our results show that the disk is likely to be strongly unstable to warping; the outcome of 
this instability is presumably the settling of the disk into the principal plane normal to the minor 
axis 9 . It is possible that some warped galaxies are a manifestation of this instability. 

Rotating halos can also excite warps. Warps in vertically stiff halos generally have retrograde 
pattern speeds (eq. 17) and positive energy, so that excitation requires that the halo and the 
disk rotate in opposite directions; as we have mentioned, cosmological simulations show that the 
angular momentum vectors of halos and their embedded disks are generally similar (although polar 
rings [Casertano et ai. 1991] and counter-rotating disks [Rubin et ai. 1992] offer counter-examples). 
Warps can also be excited if the halo orbits are predominantly circular, but simulations of galaxy 
formation suggest that halo orbits are radial rather than circular. Thus excitation by either of 
these mechanisms is unlikely to be the major cause of warps either. 

A general difficulty with explaining warps by excitation from dynamical friction is that the 
excitation timescales are short compared to a Hubble time, and there is no obvious nonlinear 
process to limit the amplitude of the warps at their observed values. 

In the majority of galaxies, for which dynamical friction damps rather than excites warps, 
the damping time is short compared to the Hubble time. In the plausible case that halos have 
predominantly radial orbits, the damping time is often shorter than even the precession time 



These results bring into sharp focus the problem of exciting warps. Ostriker & Binney 
(1989) have argued that normal galaxy formation processes will not excite the Laplace mode; our 
calculations show that even if the mode were excited during galaxy formation, it would by now 
have been damped away. Thus we require an excitation mechanism that operates throughout the 
life of the galaxy. 

One plausible mechanism is twisting of the halo by cosmic infall (Ostriker & Binney 
1989). Both analytic arguments and N-body simulations show that the direction of the angular 
momentum vector of a halo — and hence the orientation of its principal planes — is reoriented by 



A corollary is that massive polar rings should be unstable. 



2Tr/\n p \ (Figure 2). 
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infalling material with each doubling of the halo age (see Binney 1992 for a review). This steadily 
changing orientation is communicated to the inner halo and disk by gravitational torques; the 
warp arises because it is necessary to transport off-axis angular momentum between the inner and 
outer disk. In this model the warp amplitude is determined by the present rate of twisting rather 
than the past history of excitation and damping. 

Our results can be applied to warps with other azimuthal wavenumbers; although we have 
done no calculations we can make some comments. Sparke (1994) has argued that axisymmetric 
(to = 0) 'bowl-shaped' modes may be present in spiral galaxies and polar rings. The analysis of 
Section 2.2 shows that all to = warps have positive energy; hence they will be damped in an 
iddf halo. Symmetry arguments show that halo rotation will not tend to excite axisymmetric 
modes either. We expect the damping from dynamical friction to be rapid, since the frequency of 
the to = mode is usually much larger than the pattern speed of the Lagrange mode. Thus it 
is unlikely that to = modes are excited by dynamical friction; and if they are excited by other 
mechanisms, they will be rapidly damped. For higher azimuthal wavenumbers, to > 2, there is a 
continuum of bending modes in realistic disks, so such modes will disperse even without damping 
from the halo. 

We close this section with a summary of the approximations we have made: 

1. We have used linear perturbation theory (in the warp amplitude), which is inadequate: 
observed warps may have amplitudes ^ 10°, which represents only the difference in 
inclination between the inner disk and the outer edge. In perturbation theory the tilted 
disk tends to align with the equatorial plane of the halo; in reality a tilted disk with mass 
comparable to the inner halo has much more angular momentum than the halo, so that the 
inner halo tends to align with the disk rather than the disk with the halo (then over a longer 
timescale the outer halo continues to exert significant frictional force and eventually brings 
all three components to a common alignment). 

2. In calculating the orbits of halo stars, we have assumed that the unperturbed potential of 
the halo plus disk is spherical, whereas in fact the halo potential is probably triaxial and 
the disk provides an additional flattened component. With a more realistic potential the 
friction would probably be increased, since the halo would generally be flattened towards 
the disk plane so its interaction with the warp would be stronger. An inconsistency in our 
calculations is that we assume the halo is axisymmetric but non-spherical when computing 
the warp energy and pattern speed. 

3. We have approximated the halo potential as that of a singular isothermal sphere. This 
approximation leads us to underestimate the damping/excitation for warps with small 
pattern speeds, since the absence of a core implies that fewer halo stars are close to the 
(/ l5 / 2 ) = (1, —2) resonance that dominates the friction at small pattern speeds. 

4. We have neglected the self-gravity of the halo response in order to simplify the friction 
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calculations. Including the self-gravity of the halo response should generally increase the 
dynamical friction (for satellites, the inclusion of the response self-gravity decreases the 
friction [Hernquist & Weinberg 1989], but this reflects mainly the difficulty of treating the 
motion of the center of mass consistently without self-gravity), and therefore strengthen our 
conclusion that the damping times are short. 

5. We have treated the disk as razor-thin, with no internal structure. The dispersion relation 
for bending or corrugation waves in stellar systems with slab geometry (uniform in x and 
y, self-gravitating, and symmetric about the z = plane) has been investigated by Araki 
(1985), Weinberg (1991), and Toomre (1966, 1994). The bending wave is damped by 
near-resonant stars; if the horizontal variation of the bending wave is oc exp[i(k • x — ut)], 
the resonant stars are those with k • v = nv(E z ) + u, where n is an integer and v{E z ) is the 
vertical oscillation frequency as a function of vertical energy (this approximation neglects 
the curvature of the epicyclic orbits in a real galaxy disk, but this should not invalidate 
our analysis since the epicycle frequency is small compared to v). If the typical horizontal 
velocity dispersion is <7, the criterion for significant damping — plenty of resonant stars — is 
ka /v(Qi) > 1 + |w|/z/(0). The term |cj|/z/(0) is usually negligible for Laplace modes, which 
have small pattern speeds. Taking A = 2ir/k = 10 kpc, a = 40km s _1 , z/(0) = 3 X 10~ 15 s _1 
(corresponding to a local density of O.15M pc -3 ), we find ka /v{Gi) = 0.3; which is probably 
too small for significant damping in a typical warp. 

Many of the concerns raised by these assumptions can be addressed by N-body simulations, 
which complement the analytic approach described here. Dubinski & Kuijken (1994) have recently 
examined the evolution of tilted disks (both rigid and N-body) embedded in N-body halos. They 
find that the disk rapidly aligns with the inner halo, and the damping times that they observe are 
even shorter than the ones we calculate (probably because the disk contains most of the angular 
momentum, so that the halo aligns with the disk rather than vice versa). Our results suggest 
other worthwhile N-body experiments. It would be interesting to confirm whether warps can be 
excited in rotating or vertically soft halos, to test our prediction that the damping is strongest in 
halos with predominantly radial orbits, and to investigate the fate of massive polar rings. 

Our conclusions are summarized briefly in the abstract. 

We acknowledge helpful discussions with James Binney, Omer Blaes, Ray Carlberg, John 
Dubinski, Peter Goldreich, Linda Sparke, and Alar Toomre. This research was supported by 
NSERC. 

A The Dimensionless Energy Loss Rate 

We wish to evaluate the rate of energy transfer to the disk given by equation (62). 
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An orbit in a spherical potential can be described by the canonical actions 
(Ii, 12,13) = (I r ,J,J z ), and their conjugate angles (wi, W2, W3); I r is the radial action, and 
J and J z are the magnitude of the total angular momentum and its component along e z which 
we take to be normal to the unperturbed disk. The conjugate angles vary from to 2ir, with 
frequencies given by Hamilton's equations (51). 

Following Tremaine & Weinberg (1984) and Weinberg (1986), we expand the perturbing 
potential from the warped disk in both spherical harmonics and an action- angle Fourier series (cf. 
eq. 55) 10 , 

CO / 

U w (r,t) = EE U ln (r)Y ln (0,O)e m ^- Qpt) 

1=1 n =-l 

{00 00 I 

E E vim&v^-^} (ai) 
/ 3 =0/i,/ 2 = -oo J 

where I3 = \n\, u\ = h^lp, U* n = ( — l) n Ui t - n , and the radial functions JJ\ n are derived in Appendix 
B. The amplitudes in the two lines of equation (Al) are related by 

00 / 9 \ 

y hhh (hj2j3) = E [y^-j v 'hh(m!i 2h (hj2), (A2) 

where 

Vu 2l M = rU(mh(h,^ h ~ h , (A3) 

WhAl^h) = ~ I* d Wl cos[l lWl - l 2 (il> - w 2 )]Uu 3 (r); (A4) 
vr jo 

here ip is the angle in the orbital plane measured from the ascending node, and explicit formulae 
for the angles w\ and w 2 are given in Tremaine & Weinberg (1984) and below (eqs. A10 and All). 
Here cos fi = J z /J is the cosine of the orbital inclination, and rj ^ (/?) are rotation matrices which 
satisfy the orthogonality condition (cf. Gottfried 1979, section 34.5) 

q d{3 S m{3r l mn (]3)ri n (]3) = ^^S w . (A5) 

It can be shown that = 

For the logarithmic potential given in equation (68) it is convenient to define a dimensionless 
radius and angular momentum (Weinberg 1986) 



10 When n = the time-dependence must be written as e ,wt rather than e p r simplicity, we shall not 

explicitly include this special case in our formulae. 
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where J C (E) = J s e^ E l v ^ x / 2 ) is the maximum angular momentum possible for a given energy (a 
circular orbit) and J s = r s v c . In these variables the orbital frequencies can be written in the form 



ni(E,J) = e- E ^Ui(K), 
where the new frequencies are given by V{ = (v c /r s )vi, 



z/i(k) 



1 



dr 



2vr J (-21nf - K 2 /ef 2 y/ 2 
v\(k)k f dr 



(A7) 



(A8) 



(A9) 



27re 1 /2 J P(-2hir- K 2 /epy/ 2 ' 

and z/3 = 0. The integrations are taken over a complete orbit with turning points determined by 
the roots of the denominators. 11 



The angles are defined as a function of the dimensionless radius by 

wi(r) = Mk) 



v 2 {n) k 
w 2 - V = — — w i 



p (-21nf - k 2 /^ 2 ) 1 / 2 ' 
dr 



vi(k) e 1 / 2 Jf p f 2 { — 21nr — K 2 /ef 2 ) 1 / 2 



(A10) 
(All) 



We divide the directional derivative of the df in equation (70) into two parts, even and odd 
in cos /3, 



1- 



where 



and 



dF_ 
01 



F X (E,J) = F J 2 ^-%- E ^ [~ + y[2( 7 - l)/ 2 + «/ 3 ]} , 



F 2 (i?,«7) = aFo«7 2(7 - 1) e- S/ff2 



(27-3), i-n 

h s~ 



(A13) 



(A14) 



J <T< 

Here we have used ft = dE/dl. Using dl 2 dl^ = JdJd(cos /3), the integral over (3 in equation (62) 
due to the first term, F\(E,J), is proportional to 



d/Jsin/^I)! 2 ^ 



^(2/ + l)(l + ^ 3 o) : 



■\Yi h (^,o)\ 2 \wll 2h (ii,h)\ 2 , (Ais) 



1 Weinberg sometimes uses dimensionless units r s = a = 1; in these units v c = 2 



1/2 
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which follows from equations (A2) and (A5). The integral arising from the second term, 
F 2 (E, J) cos /3, is proportional to 



[ ^sin/Jcos^I)! 2 = ^ J2 ^J/!,^)^^^/!,^)^^^^)^^^^) 



X f #isin/3cos/3r' , (P)r l ' h (fi). (A16) 



The last integral vanishes unless / = /' ± 1 (cf. Gottfried 1979); we show in Appendix B that the 
potential contains only terms odd in / + I3 so this integral will vanish. Consequently, using (A15) 
in (62) we find 

E W = Y. (2/+ l)(l + ^ 30 )2 l y " 2 (2^ ? 0)| 2 J dhJdJlWhfF^E, J)S(\ ■ n - ui). (A17) 
The resonance constraint can be written as a delta function in energy, 



(A19) 



where 

E l (K) = v 2 c ln ^ 

L ^1 

is the resonant energy associated with a particular k and is the unit step function. We replace 
dl\ by dE/Cli(E,J); then using equations (71), (A6), (A7) and (A13) we find 

JdJF 1 (E, J) 2F J^ k 2 ^" 1 ) 



where 



f(K)dK (A20) 



/(«) = 7 «1 • £/(«) - e 1 / 2 [( 7 - l)/ 2 + \ah}. (A21) 

and we have replaced the frequencies V{ by their dimensionless counterparts V{. After integration 
over energy, the integral in equation (A17) becomes 

wi f dI 1 JdJ\Wjl 2h \ 2 F 1 (E,J)6(l-n-u l ) = 

where the dimensionless transform is (cf. eq. A4) 

W l l f 2h {K^ l ) = - J d^co S [hw 1 -l 2 (iP-w 2 )}Uu 3 l K -^-), (A23) 
wi = iW^c and C/„ 3 (r) = (GM d / R d )U lh (r / R d ). 
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Finally, inserting this result in equation (A17) and using equation (71) we obtain the form we 
use for numerical computations, 

E w = ^§^m p ;~t,a), (A24) 

where 

oo oo oo , 

5(fi p ; 7 ,a)=E £ ^-r^-l^aivr^)! 2 ^^^^), (A25) 
; 3 =o h,l 2 =-oo l=i ^ > 

and 

Sf, (12 p ; 7 ,a) = - " g ( : p l 1 , / rfK^/Wl^.^^j O F . (A26) 

" 2 ' 3V p ' " 7 (l + ^ 3 o) 2 e"T( 7 ) io ^i(k) V 71 " 2 ' 3V Pn \l 3 nj V 7 

When 7 = 1 and a = 0, the df is that of a singular isothermal sphere. In this case, we recover 
the result of Weinberg (1986; eq. 40), after correcting two minor typographical errors (the factor 
21 + 2 in Weinberg's paper should be 21 + 1; and the sum over I3 should run from to /, not — / to 
/) and converting from torque to energy loss. 



Note the symmetry 

^(^W,") = ^;\ i3 (-fi p ; 7 ,-a), (A27) 

which implies that 

£(fi p ; 7 ,a) = £(-fi p ; 7 ,-a). (A28) 

Thus if the halo is non-rotating (a = 0) and I3 ^ 0, E w is independent of the sign of the pattern 
speed as one would expect from symmetry. 

The term in equation (A21) proportional to the halo spin a contributes energy to the bending 
wave if the pattern speed is in the same direction as the halo rotation, and removes energy if the 
pattern speed and halo rotation are opposite. 



B Gravitational Potential of a Disk Bending Mode 



The potential of a thin, warped disk of surface density fJ,(R) is 



U d (r,t) 



G / <ZR'- 



r- r'(R',i) 



where r' = R' + Z(^R' ,t)e z . Subtracting off the potential of the unperturbed (flat) disk, and 
expanding to first order in the small displacement Z, we obtain the potential due to the warp 



U w (r,t) 



-G [ dR'^R^ZiR^t)^—^ 
J oz' r — r 



:b2) 



z'=0 



- 34 - 



In spherical coordinates 



r — r 



In ^> 



so that 



d 1 



E 



47r r\ 



dz' |r — r' 

The values of the spherical harmonics are 

QYln 



1 dY? n n 



R' dcos0 



vr,0)e" 



-incf)' 



(B3) 



(B4) 



I, o-r-i^- 1 )/ 2 f^Y^E ^ (*+*)" 

dcos# V2 ' J 1 ' V 4vr 7 y(/ + n )!(/_ n _ l)!! 

where the double factorial terms are 

(2nV 

(2n-l)!!=^ — J i, (2n)\\ = 2 n n\. 

Terms with / + n even vanish. 

We now substitute equation (7) for Z(R, /). Then we may write 

U(r,t) = 52U, n (r)Y, n (0,Oy( n *-»» t \ 



n + / odd 



(B5) 



(B6) 



l.n 



where 

Uin(r) 
and uj m = u, LO— n 



^ 2 G dYf m 

21 + ldcos0 [2 " , " J Jo ^ y+i 



(ivr, 0) / dR' n{R')-fc[h{R')6 mn + h*(R')S 
> 



m, — n\ •} 



(B7) 



-or 



In this paper we mostly consider disturbances with m = 1. Since n = +m and n + / must be 
odd, only terms with even / will be non-zero. We assume that the disk surface density and warp 
shape are given by equations (77) and (79). Then equation (B8) can be evaluated as 



tt ( ^ GMd 27r dY *±i n ^ r , f , ^ xl < -x 



(B9) 



where 



dxf(x) 



4 +1 



C 



RA l+1 l(n + I + 1,(1 + X)rjR d ) 
r J (l + X) n + l + 1 

fjTS 1 T(n-l,(l+X)r/R d y 

\Rd) (i + \y- 1 



In this equation the gamma functions are defined by 



~/(p,z) = / x p 1 e x dx, T(p,z) 
Jo 



(BIO) 



(Bll) 
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note that 

1 (p,z) + T(p,z) = T(p) = (p- 1)1 (B12) 

If p is a positive integer then 

T(p,z) = (p-l)\e- z Y j Z -^. (B13) 
k=o K - 

For negative integers, — n < 0, 

r(-n,z)=^£„ +1 (z), (B14) 
where E n (z) is an exponential integral. 

C The Contribution of Near- Radial Orbits to Dynamical Friction 

In this Appendix we elucidate the reason for the apparent divergence in the damping/excitation 
rate in halos with predominantly radial orbits. 

Suppose that a torque er(/) is exerted on the star in a nearly radial orbit. The resulting 
change in the angular momentum vector may be written 

AJ = eA, (CI) 

where A = J r(t)dt. The change in J = |J| is then 

~ A 2 (J -A) 2 ' 



J A 2 

AJ = e— + e 2 



O(S). (C2) 



2J 2J 3 

As J 0, in general A approaches a fixed, non-zero value. Thus as J 

e 2 A 2 sin 2 -0 . 
AJ -> eA cos ^ + — — -, (C3) 

where ip is the angle between A and J (plus terms that are 0(e 3 ), which we neglect). The 0(e) 
term can be dropped, since it averages to zero in a spherically symmetric distribution of stars, and 
in any case does not diverge as J —> 0. In a spherical distribution, the average of sin 2 ip is ^, so 
the average change is (AJ) = |e 2 A 2 /J. The corresponding energy change is 

(AE) = f^(AJ) = fi 2 (AJ), (C4) 

which also diverges as J -1 , since &2 is generally constant and non-zero as J — ► 0. 

The number of stars in a small interval of energy and angular momentum is 
dN = 16ir 3 F(E, J)dE JdJ/il-y. Now assume that the df has the form (70), that is 
F(E,J) = J 2 ( 7_1 ) g(E). Since fii is constant as J — ► 0, the total energy change of the 
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stars of a given energy on near- radial orbits is proportional to / J 21 ' 1 (AJ)dJ = \e 2 A 2 J J 2 T 2 dJ, 
which diverges if 7 < \. 

Thus we have explained the apparent divergence in dynamical friction for dfs with 7 < \. 
Our arguments also show that the divergence is an artifact of the expansion procedure: the exact 
change in angular momentum is 

A/exact = (J 2 + 2eJ • A + e 2 A 2 ) 1/2 - J, (C5) 

and the total energy or angular momentum change from near-radial orbits (J < «7 max , say) is 
proportional to 

/ J^-^AJexact)^- (C6) 

JO 

When 7 > \ this integral is 0(e 2 ) (after averaging over orientations assuming (cost/') = 0); this 
is the usual case in which the frictional energy or angular momentum change is second-order in 
the perturbation strength e. When 7 < \ the integral is 0(e 27+1 ); thus the frictional force is of a 
lower order in the perturbation strength than e 2 , but is not divergent. 
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